Exploring county-level spatio-temporal patterns in opioid overdose related emergency department visits

Opioid overdoses within the United States continue to rise and have been negatively impacting the social and economic status of the country. In order to effectively allocate resources and identify policy solutions to reduce the number of overdoses, it is important to understand the geographical differences in opioid overdose rates and their causes. In this study, we utilized data on emergency department opioid overdose (EDOOD) visits to explore the county-level spatio-temporal distribution of opioid overdose rates within the state of Virginia and their association with aggregate socio-ecological factors. The analyses were performed using a combination of techniques including Moran’s I and multilevel modeling. Using data from 2016–2021, we found that Virginia counties had notable differences in their EDOOD visit rates with significant neighborhood-level associations: many counties in the southwestern region were consistently identified as the hotspots (areas with a higher concentration of EDOOD visits) whereas many counties in the northern region were consistently identified as the coldspots (areas with a lower concentration of EDOOD visits). In most Virginia counties, EDOOD visit rates declined from 2017 to 2018. In more recent years (since 2019), the visit rates showed an increasing trend. The multilevel modeling revealed that the change in clinical care factors (i.e., access to care and quality of care) and socio-economic factors (i.e., levels of education, employment, income, family and social support, and community safety) were significantly associated with the change in the EDOOD visit rates. The findings from this study have the potential to assist policymakers in proper resource planning thereby improving health outcomes.


Introduction
The continued rise in drug overdoses involving opioids has significantly impacted the social and economic fabric of communities in the United States [1]. For instance, in 2017, the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 economic cost of opioid deaths, criminal justice involvement, treatment, and lost wages surpassed $1.02 trillion [2]. Between April 2020 and April 2021, over 100,000 deaths involved opioids-an increase of 28.5% over the previous 12 months [3]. An effective response to the opioid crisis requires policymakers to understand which communities are the most affected and how the resources have been allocated in those communities [4].
Recent research has identified many socio-ecological factors as the risk factors of opioid overdose [5,6]. These risk factors are inextricably linked to personal and environmental factors and systems that hinder rather than support individuals [6]. For instance, among individuals receiving healthcare services at a free clinic, prescription opioid misuse was more likely among patients who were employed and less likely among those with post-high school education [7]. Uninsured individuals were significantly more likely than insured individuals to be high-risk drug users [8]. Access to prescription opioids is another prominent risk factor of unintentional opioid overdose deaths [9,10], which has increased in rural areas with a greater need for medical services [11,12].
Opioid overdose rates within the United States have been found to vary across different geographical regions. A socio-ecological framework posits that the characteristics of the community in which individuals live significantly influence their health behaviors [13]. For instance, counties across the United States with high economic distress, high opioid prescription rates, and a lack of opioid treatment program providers have higher opioid overdose mortality rates [14,15]. Statewide health disparities, including lower socioeconomic status and access to health care, can differ at the county level and are more predominant in rural areas [16]. Opioid overdose patterns are consistent with counties that may lack the resources necessary to prevent overdose [17]. In sum, the intersection of where and how people live is a significant factor in health outcomes and requires public health to work with urban planning to create supportive and healthy environments to reduce the risks of opioid overdose [18].
Besides geographical differences, there have been significant variations in the temporal trends of opioid overdose rates. Beginning around the year 2000, opioid overdose rates in the United States have been increasing over time [19]. In recent years (since the early 2010s), much of this growth has been attributed to synthetic opioids such as fentanyl, which increased more than 50% from 2019 to 2020 [20]. During the same time period, prescription opioidrelated overdose showed the first increase in years (10.6%) whereas heroin-related overdose showed a downward trend (down 3.6%), similar to the recent prior years [20]. Not all regions within the United States have the same characteristics and the national pattern may not well reflect the local growth trajectories.
Prior research has implemented different spatiotemporal analysis techniques to identify the local geographical differences in opioid overdose rates over time. For instance, Hernandez et al. [21] examined prescription opioid death rates in Ohio from 2010-2017 and identified 12 hotspots along with three significant changing trends of opioid overdose using temporal trend analysis. Marotta et al. [17] examined cumulative opioid overdose deaths in New York State using data from 2013-2015 and identified geographical hotspots of overdose death rates for different types of opioids. Sauer et.al. [22] used spatio-temporal bayesian modeling and exploratory spatial analysis to evaluate risk factors related to drug-involved emergency department visits in the greater Baltimore metropolitan from 2016-2019. These studies and a few others [23][24][25] have demonstrated that spatio-temporal techniques utilizing local population-level data can provide a profile of opioid overdose risk.
In this study, we performed a county-level spatio-temporal assessment of opioid overdose rates and their association with different socio-ecological factors for the state of Virginia. We focus on Virginia for two main reasons. First, its growing rates of fatal opioid overdose [26] represent the growing rates of opioid overdose across the United States. Secondly, the Virginia Department of Health collects monthly emergency department opioid overdose (EDOOD) visit data as part of syndromic surveillance to measure health trends [27]. This publicly available dataset can serve as a timely indicator of opioid overdose trends. Different from prior studies, we used the EDOOD visits as a proxy for total overdoses to understand the spatio-temporal dynamics of opioid overdose rates and potential socio-ecological risk and protective factors. Emergency departments (EDs) are the primary treatment venue for patients with overdoses [28]. In recent years, urgent care centers' utilization to treat opioid overdoses has also significantly increased. Between 2007-2016 documented claims for urgent care centers increased by 1,725 percent compared to a 229 percent increase for emergency department claims with 'injury, poisoning, and consequences of external causes' [29]. However, the current study only focuses on the visits obtained from hospital-based and free-standing EDs. EDOOD visit rates increased by 28.5 percent across the United States in 2020, compared to 2018 and 2019 [30]. Understanding the spatio-temporal trends of EDOOD visit rates can help identify targets for policy change and timely resource allocation to mitigate the opioid crisis.
To summarize, this study utilized a comprehensive, three-pronged approach to understanding the opioid overdose trends in Virginia. The goals of the study were to 1) identify spatio-temporal variations of EDOOD visit rates from 2016-2021 among Virginia counties, 2) assess how counties cluster together based on their EDOOD visit rates, and 3) identify socioecological factors that are associated with the change in EDOOD visit rates over time. Moran's I [31], Local Indicators of Spatial Association (LISA) [32], Dynamic Time Warping (DTW) [33], and multilevel modeling [34] were implemented for the spatio-temporal analysis. To our knowledge, this is the first study to combine techniques from statistics, data mining, and geographic information systems (GIS) to examine how a county performs in terms of EDOOD visits. Although the study focuses on Virginia, the study methods can be extended to other geographic locations with similar data. The source code can be made available upon request.

Study area
We analyzed the EDOOD visit rates and the associated socio-ecological factors across the different counties within the state of Virginia, United States. The state of Virginia consists of 95 counties and 38 independent cities that are considered county-equivalent for census purposes. The analysis was performed for those 133 unique geographic regions.

Measures
Emergency department opioid overdose visits. The EDOOD visits dataset was obtained from the Virginia Department of Health (VDH) [35] and is based on syndromic surveillance reported by hospitals and free-standing EDs in Virginia. It consists of count and rate statistics (monthly and annual) of ED visits for unintentional opioid overdose (fatal and nonfatal) among Virginia residents aggregated by different geographical units. This dataset excludes heroin from all other types of opioid overdoses. There is separate data that reports on heroinrelated ED visits. In Virginia, only a small percentage of the total overdoses are driven by heroin in recent years [35]. The numbers are not significant enough to analyze the spatio-temporal variations of opioid overdose rates. Thus, this study did not include data on heroin-related ED visits in its analyses.
Although the EDOOD data spans from 2015 to 2022, we only examined data from 2016 to 2021 given that complete data (i.e., previous 12-month average visits, annual summary) was not available for 2015 and 2022. The outcome variable was defined as the rate of EDOOD visits per 100,000 population.
Geography assignment. The county location in the EDOOD visit dataset is assigned based on the patient's self-reported residential zip code. A single zip code may belong to multiple cities and counties. In that case, the visit is assigned to the county/city where the majority of the population resides. Additionally, some Virginia cities and counties are aggregated (e.g., Alleghany County and Covington City) due to overlapping zip codes (see S1 File).
EDOOD visit definition. ICD9 (e.g., 965.00, 965.01), ICD10 (e.g., T40.0X1A, T40.601A), or SNOMED (e.g., 295165009, 242253008) codes representing opioid overdose were used to identify overdose incidents. Similarly, the mention of terms like Narcan or naloxone in the chief complaint or discharge diagnosis was also used to identify overdoses. This includes unintentional overdose by opioids or unspecified substances (excluding heroin). The definitions mentioned here are merely provided as examples. For complete information on the inclusion and exclusion criteria, please refer to the "Unintentional overdose by opioid or unspecified substance (excluding heroin)" section in [36].
Converting visit counts to rates. Whenever information was only available in the form of counts, we used the population data from American Community Survey (ACS) to calculate the rates (visits per 100K population) [37]. We followed the same guidelines as specified in . We extracted data from 2016 to 2021 to be consistent with the EDOOD visit dataset. This publicly available dataset aggregates data from the Centers for Disease Control and Prevention (CDC) as well as other sources (e.g., U.S. Census Bureau, Behavioral Risk Factor Surveillance System) to provide yearly county-level rankings based on different attributes related to health outcomes. Based on prior studies on the association between socio-ecological factors and opioid overdoses [5-9, 14, 15], we selected four variables: health behaviors, clinical care, social and economic factors, and physical environment as the most appropriate for our study (see Fig 1). The goal was to examine if these socio-ecological factors influence the spatiotemporal trends of EDOOD visit rates. Other sources have also aggregated data on socio-ecological factors (e.g., Opioid Environment Policy Scan, Virginia Department of Health, etc.). However, the CHR&R dataset is an adequate proxy for socio-

PLOS ONE
County-level spatio-temporal patterns of opioid overdoses ecological factors for our analysis because it encompasses a wide range of socio-ecological inputs within its four variables: Clinical care. Includes inputs about access to care and quality of care, such as uninsured rates, and access to primary care physicians, dentists, and mental health providers.
Social and economic factors. Includes inputs from six unique data sources measuring unemployment, children in poverty, income inequality, single-parent households, social associations, violent crime, and injury deaths.
Physical environment. Includes inputs from four unique data sets measuring air pollution, alcohol drinking violations, severe housing problems, driving alone to work, and long commute-driving alone.
Health behaviors. Includes inputs about tobacco and alcohol use, diet and exercise, and sexual activity from seven unique data sources. These inputs further include physical inactivity, excessive alcohol use and impaired driving deaths, sexually transmitted diseases, and teen pregnancy.
The rankings in the CHR&R dataset were calculated using standardized z-scores from several data sources [39]. The county with the lowest z-score received a rank of 1, which indicates the highest quality of socio-ecological factors (e.g., low tobacco use, better access to care, better education). The county with the highest z-score received a ranking of 133, which indicates the lowest quality of socio-ecological factors (e.g., high tobacco use, poor access to care, poor education). The z-scores were averaged for the counties (or cities) that were combined in the EDOOD visit dataset. The CHR&R dataset has been validated in other studies [40,41].
Neighborhood adjacency. Data on neighborhood adjacency for Virginia counties was obtained from the US Census Bureau [42]. This data lists each county along with its adjoining neighbors including counties that are not in Virginia but adjacent to Virginia counties. For this study, we only considered the neighboring counties that are a part of Virginia. This data was utilized for our spatial analysis.

Data analysis
Spatial analysis. To identify any spatial variations in the opioid overdose rates across Virginia counties, we calculated the spatial autocorrelation using data on neighborhood adjacency and EDOOD visits for the years 2016-2021. Average monthly EDOOD visit rates were used to summarize the yearly EDOOD trends. Spatial autocorrelation is the phenomenon where the presence of some quantity in an area makes its presence in neighboring areas more or less likely [43]. Positive autocorrelation, which is more common in practice, is the tendency for areas that are close together to have similar values. In contrast, negative autocorrelation is the tendency for areas that are close together to have different values.
Global spatial autocorrelation. Global autocorrelation measures the overall association within the data. In this study, it measures the similarity between the neighboring counties in terms of EDOOD visit rates. We calculated the Moran's I index, a common measure of global autocorrelation [31]. Then we performed a permutation test to assess the significance of Moran's I index analysis. The values of Moran's I range from +1 (strong positive spatial autocorrelation) to 0 (randomness) to -1 (strong negative pattern). Moran's I value of 0.7, for instance, indicates that the spatial pattern across counties is homogeneous meaning that the neighboring counties have very similar visit rates.
Local spatial autocorrelation. To study the contribution of each county to the global Moran's I index and identify local hotspots (clusters of high EDOOD visit rates) and coldspots (clusters of low EDOOD visit rates), we calculated Local Indicators of Spatial Association (LISA) [32] for each county. These auto-correlation indices were used to divide the counties into four distinct groups: • High-high: Counties with high visit rates with neighboring counties that also have high visit rates (also known as hotspots) • Low-low: Counties with low visit rates with neighboring counties that also have low visit rates (also known as coldspots) • High-low: Counties with high visit rates but surrounded by counties that have low visit rates • Low-high: Counties with low visit rates but surrounded by counties with high visit rates The classification of counties into regions of low or high visit rates was done based on whether they had rates less than or greater than (or equal to) the mean value of visit rates across the state of Virginia. A permutation test was used to identify non-significant associations within neighbors (counties) thereby highlighting the significant hotspots and coldspots of EDOOD visit rates. Both local and global autocorrelation analyses were performed in python using the PySAL package [44].

Temporal analysis
Temporal analysis refers to the study of an outcome over time. We used two different methods -clustering and multi-level modeling-to analyze the EDOOD visit rates over time and identify their association with different socio-ecological risk factors.
Clustering. To identify similarities between the temporal trends of visit rates in different counties, we used dynamic time warping (DTW) [33]. DTW is a data mining technique used to compute the similarity between multiple time series (e.g., opioid overdose rates over time) and cluster them together based on their shapes and magnitudes. We utilized the moving average of the monthly visit rates (for a smoother curve) to perform the clustering. After we obtained clusters of counties, we mapped the clusters using choropleth maps for easy visualization. Analyses were performed in Python using the PySAL package [44] and tslearn package [45].
Multilevel modeling. To identify how the changes in EDOOD visit rates relate to the changes in different socio-ecological factors, we performed multilevel modeling [34] with the visit rates as the outcome variable and four different time-varying aggregated variables from the CHR&R dataset as our predictors. Multilevel modeling is advantageous because it accounts for correlations across time within individual counties. Moreover, multilevel models handle missing data in visit rates for any county and at any time point without pairwise deletion of individual counties. Multilevel modeling was performed in R using the lme4 package [46]. A forward-stepping procedure was used to create the final model [47]. First, an unconditional means model (i.e., baseline model) was created. From this model, an intraclass correlation coefficient (ICC)-representing the proportion of variance explained within counties-was computed. Next, conditional growth models were created to examine the linear effect of time on EDOOD visit rates, with time modeled as a fixed and random slope in separate models. The model (i.e., either fixed or random slope) with a better fit compared to the unconditional growth model was used moving forward. Finally, a conditional random growth model was created to examine the linear effects of the time-varying covariates (i.e., socio-ecological factors) on the visit rates: Level 1 :

PLOS ONE
County-level spatio-temporal patterns of opioid overdoses Level 2 : The Level-1 equation models the within-county variance based on EDOOD visit rates. Thus, for county i at time t, the expected outcome, Y, is equal to the intercept, π 0i , plus an effect for the slope, π 1i , plus an effect for Health Behaviors, π 2i , plus an effect for Clinical Care, π 3i , plus an effect for Social and Economic Factors, π 4i , plus an effect for Physical Environment, π 5i , plus error, e it . The Level-2 equations state that the intercept and year were fitted using random effects whereas the socio-ecological predictors were modeled using fixed effects. We used a fixed slope for the socio-ecological factors because specifying many random coefficients overfits the model producing misleading results [48].
Lastly, we ran ANOVA-like table with tests of random effects for each model using the ImerTest package. Each model was compared to the preceding model using these ANOVA tests.

Spatial analysis
Global auto-correlation. The results from Moran's I index indicate that there is some similarity between counties and their neighbors with respect to their EDOOD visit rates in most of the years. As shown in Table 1, the values of indices are greater than 0 for all the years indicating a positive global spatial autocorrelation. However, the similarity scores and the corresponding p-values vary across different years. The strongest association seems to be present in the year 2018 (Moran's I = 0.25). This indicates that many neighboring counties had similar EDOOD visit rates during that time. On the other hand, the Morans'I value for 2016 is almost close to 0, meaning that the EDOOD visit rates were not significantly similar across the neighboring counties.
Local auto-correlation. Fig 2 represents the plots of EDOOD visit rates alongside the classification returned by LISA for the 3 years (2017, 2018, and 2021) with the most significant neighborhood association (based on Moran's I values). To show them in the maps, we assigned the same value/color coding to counties that were combined together in the analyses (although they were only counted once in analyses). For instance, the combined rate for Grayson and Galax for the year 2016 was 21.7, but they are both represented separately as a value of 21.7 on the map.

PLOS ONE
County-level spatio-temporal patterns of opioid overdoses EDOOD visit rates. As shown in Fig 2A, The EDOOD visit rates differed within the state and across time. However, some of the highest rates of visits can be seen across the southwestern region (e.g., Galax & Grayson, Smyth, Martinsville & Henry) and also in the northwestern region (e.g., Orange, Louisa, Culpeper). In contrast, most counties in the northern region (e.g., Fairfax & Falls Church, Loudon) and some in the eastern region (e.g., Accomack) had relatively lower rates.
LISA values. The LISA values were calculated for every county in Virginia for the years 2016-2021. However, as aforementioned, Fig 2B only presents results for the years 2017, 2018, and 2021. The four distinct subgroups returned by LISA are shown in choropleth maps (in Fig  2B) in distinct colors: red (high-high), blue (low-low), light blue (low-high), pink (high-low). The significant clusters returned by the permutation test are highlighted in yellow. These highlighted counties had significantly higher or lower concentrations of EDOOD visit rates based on their color coding (red: hotspots, blue: coldspots). As expected, the hotspots were mostly concentrated around the southwestern region and the northwestern region. The cold spots were scattered across multiple regions with the most consistent one being in the northern region. There were some counties (pink clusters and blue clusters) that had significantly different visit rates than their neighboring counties. For instance, counties like Floyd, Scott, and Caroll in southwestern Virginia had lower EDOOD visit rates even when most of their neighboring counties had higher rates. Similarly, as also verified by Moran's I, the neighborhood similarity seems to be the most prominent in the year 2018, where there are close clusters of high and low EDOOD visit rates. The locations of the LISA subgroups seem to be changing over time.

Temporal analysis
Clustering. Five distinct groups of counties (clusters) were returned by the DTW clustering algorithm. The clusters with plots of their temporal trends and their corresponding geospatial mapping are provided in Fig 3. These clusters differ in their magnitude and trend over time. As depicted in Fig 3, counties in Groups A and E had similar trends over time but differed in their magnitudes. Most counties in these groups had EDOOD visit rates that slightly decreased from 2016 to 2018, which started increasing around 2019. However, counties in Group A (e.g., Buchanan, Louisa, Shenandoah) had slightly lower rates of EDOOD visits (starting range 7-16) as compared to the counties in Group E (e.g., Orange, Smyth, Wise, Dickenson), which had a starting range of 13-20. Similarly, counties in Group B (e.g., Amelia, Bland, Amherst) and Group C (e.g., Grayson & Galax, Martinsville & Henry) also had similar trends which differed in magnitudes. Overall, these counties had increasing rates over time.

PLOS ONE
County-level spatio-temporal patterns of opioid overdoses Group B had a starting range of 0-10 whereas Group C had a starting range of 0-15. Group C, however, had a steeper rising curve with an ending range of 16-38. Lastly, counties in Group D (e.g., Arlington, Fairfax & Falls Church, Alexandria, Williamsburg) did not show a consistent temporal trend and mostly had low EDOOD visit rates throughout the 6-year period. Note that for each month, we plotted the average EDOOD visit rates of the previous 12 months.
Multilevel modeling. The baseline unconditional model returned an ICC of 0.59, which indicates that 59% of the variance is attributable to differences between counties while 41% of the variance is attributable to differences within counties over time. Given that more than 5% of the variance is attributable to differences within counties over time, the use of multilevel modeling is justified. We ran ANOVA-like table tests of random effects for each model and compared each model to the preceding model. First, we found that the fixed growth model was a better fit than the unconditional model (X 2 = 34.37, df = 1, Pr(>Chisq) = 4.556 x 10 −9 , p < .001). Next, we found that the random growth model was a better fit than the fixed growth model (X 2 = 60.197, df = 2, Pr(>Chisq) = 8.481 x 10 −14 , p < 0.001). Thus, we proceeded to run a random growth model with our time-varying socio-ecological predictors. Table 2 shows the goodness of fit values for all four models. The random growth model showed that time predicted, on average, increases in EDOOD visits from 2016-2021 (see Table 2). When incorporating predictors into the model (i.e., conditional random growth model), Clinical Care and Social and Economic Factors emerged as significant time-varying predictors of the slope for EDOOD visits (see Table 2). This suggests that a 1 unit decrease in Clinical Care z-scores (i.e., higher ranking/better clinical care) increased the slope of EDOOD visit rates by 6.64. In addition, a 1 unit increase in Social and Economic Factors z-scores (i.e., lower ranking/lower quality of social and economic factors) increased the slope of visit rates by 6.74. These sociological factors are important predictors of changes in EDOOD visits.

Discussion
This study examined spatio-temporal patterns of EDOOD visits across Virginia, as well as socio-ecological factors associated with these patterns using techniques from statistics, data mining, and geographical information systems (GIS). Our spatial analysis revealed that EDOOD visit rates significantly varied across Virginia counties with clusters of overdose hotspots (primarily southwestern region) and coldspots (primarily northern region). Clustering analysis helped identify 5 distinct groups of counties based on the magnitude and the direction of change of the EDOOD visit rates over time. Although the overall trend of EDOOD visit rates differed in these groups, we observed rising trends in recent years (starting around 2019) and a slight decline in the visit rates from 2017 to 2018, in all the groups. The steepest rise in the EDOOD visit rates was seen in counties belonging to Group C (e.g., Grayson & Galax, Martinsville & Henry). Finally, the multilevel analysis revealed that the changes in the EDOOD visit rates were significantly associated with the changes in clinical care factors (i.e., access to care and quality of care) and socio-economic factors (i.e., levels of education, employment, income, family and social support, and community safety).
As aforementioned, hotspots of EDOOD visit rates were the most prominent and consistent in the southwestern part of Virginia. Southwest Virginia is a rural area where residents have higher morbidity and mortality rates often correlated to a shortage of healthcare services [49]. Residents in southwest Virginia often cannot afford annual health insurance deductibles and many medical expenses are not covered by insurance [50]. A portion of individuals reported only seeking healthcare as a last resort and many did not receive regular care from a health provider. Individuals in southwestern counties with lower access to care and quality of care, therefore, may be more at risk of opioid overdose. For example, Martinsville County had one of the highest EDOOD rates throughout the six-year period. On the other hand, many counties in northern Virginia had low EDOOD visit rates throughout the 6-year period. This region includes counties that are close to the capital and is considered to have good access to healthcare, better employment opportunities, and higher levels of education [38].
Study findings also revealed sub-groups of counties with similar EDOOD visit trends. Spatial mapping of these counties in Fig 3 indicates that many counties that belong to the same group are clustered together in space suggesting a possible association with the neighborhood characteristics. In most counties, the EDOOD visit rates decreased from 2017 to 2018. This is consistent with the national pattern and is believed to be attributed to the reduced prescribing volume of high-dose opioid pills and a sudden decline in the availability of a highly potent synthetic opioid (carfentanil) [20]. The increase in EDOOD visits from 2019 may depict how Covid-19 influenced the overall opioid overdose trends. Research conducted across six health care systems in the US [51] identified that EDOOD visit counts increased by 10.5% in 2020 compared to the counts in 2018 and 2019 despite having a 14% decrease in the overall ED visits. The same study pointed out that this rise might be attributed to the disruption of access to treatment, social support, loss of employment, social isolation, and many more during the Covid-19 period. Further, this increase in EDOOD rates from 2019 may reflect the increased prevalence and use of illicitly manufactured fentanyl, which has been the major driver of the opioid epidemic over the past few years (having surpassed prescription opioids and heroin as major causes of opioid overdose) [52]. This may also explain the sharp, increasing rates of EDOOD visits (throughout the 6-year period) in some Virginia counties (Group C in Fig 3).
Finally, our multilevel analysis found that a decrease in socio-economic factors over time was associated with increased EDOOD visit rates. This finding is consistent with studies that demonstrate associations between poor economic and social conditions and high opioid overdose mortality rates [5,6]. Conversely, the analysis revealed that counties with poor access to care and quality of care (i.e., higher clinical care z-scores) had lower EDOOD visit rates. It is possible that the inaccessibility and poor quality of clinical care might have led to individuals not being able to seek care or go to the ED. These findings suggest that addressing county-level deficits in clinical care and socioeconomic factors may help reduce opioid overdoses.
Some limitations of this study could be addressed in the future. If the necessary data is available, the study can be expanded to include a broader time frame and nationwide data. Additionally, breaking down the EDOOD visits by the type of opioids could provide a better picture of the epidemic which was not possible with the data that we used. Lastly, the spatiotemporal variations of EDOOD visits were assessed separately but they could be modeled together to identify how these interact with other outcomes and with each other.

Conclusion
Overall, there are differences between the counties in Virginia in their EDOOD visit patterns across time. These differences are significantly associated with socio-economic factors (i.e., education, employment, community safety, income, and family and social support) and clinical care (i.e., access to care and quality of care). Targeting areas that are consistently hot spots for EDOOD rates and identifying areas that vary over time is critical to address the social determinants of opioid use disorders and health care access.